function components = calculateFFT(y,Fs, NPeaks)
% Calculates the FFT and returns the main NPeaks
% [freq_1 ampl_1 freq_2 ...]

L = length(y);      % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);

components = [];
YY = Y(1:NFFT/2);
for i=1:NPeaks
    [val,pos,YY] = obtain_peaks(YY);
    peak_freq = round(f(pos));
    components = [components; peak_freq val];
end

end

function [val,pos,YY] = obtain_peaks(Y)
% Aux functions that obtains a peak and sets the 
% amplitude = 0 in the peak sorroundings
npoints = length(Y);
YY = Y;
[val,pos]=max(abs(YY));
%disp(['Peak:' num2str(val) ' ' num2str(peak_freq) ' Hz']);
peak_amplitude = 10;

p1 = pos - peak_amplitude;
p2 = pos + peak_amplitude;
if p1 < 0
    p1 = 1;
end
if p2 > npoints
    p2 = npoints;
end
for i=p1:p2
    YY(i) = 0;
end

end
